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ABSTRACT 

Evidence of jet precession in many galactic and extragalactic sources has been reported 
in the literature. Much of this evidence is based on studies of the kinematics of the jet 
knots, which depends on the correct identification of the components to determine their 
respective proper motions and position angles on the plane of the sky. Identification 
problems related to fitting procedures, as well as observations poorly sampled in time, 
may influence the follow up of the components in time, which consequently might 
contribute to a misinterpretation of the data. In order to deal with these limitations, 
we introduce a very powerful statistical tool to analyse jet precession: the cross-entropy 
method for continuous multi-extremal optimisation. Only based on the raw data of 
the jet components (right ascension and declination offsets from the core), the cross- 
entropy method searches for the precession model parameters that better represent 
the data. In this work we present a large number of tests to validate this technique, 
using synthetic precessing jets built from a given set of precession parameters. Aiming 
to recover these parameters, we applied the cross-entropy method to our precession 
model, varying exhaustively the quantities associated to the method. Our results have 
shown that even in the most challenging tests, the cross-entropy method was able 
to find the correct parameters within 1%-level. Even for a non-precessing jet, our 
optimization method could point out successfully the lack of precession. 
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1 INTRODUCTION 

Improvements in the sensitivity and angular resolution of 
astronomical instruments have allowed the study of galac- 
tic and extragalactic jets in relative detail. In our Galaxy, 
jets have been observed directly in star-forming regions 



(e.g.. IReipurth et al.ll2004 ) , plane tary nebula (e . g., Hugginsl 
120071 ') and in micro-quasars (e.g.. IStirling et al.ll2002l ). The 
Galactic Centre might also harbour a very faint jet gen- 
erated in the vicinity of a super-massive black hole (e.g., 
IKrichbaum et afll 19951) . In the case of extragalactic sources, 
jets have been detected in the cores of a ctive galaxies , some 
of them having relativistic speeds (e.g.. IZensusl 11997). Dis- 
crete components, which nature is still a matter of debate, 
are often seen receding from the unresolved core, where the 
engine respo nsible for their acceleration i s supposed to be 
located (e.g.. lJunor. Biretta fc Liviolll999l ). 
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The accumulation of observational data in a wide 
range of radio frequencies have shown that a large num- 
ber of jets are not perfectly linear, exhibiting a bent shape 
that has been interpreted as d ue to jet-ambient interac- 
tion (e.g.. iGallimore et al.|[20u4), hydro- or magnetohydro- 
dynamical instabilities (e.g-7|F crrari 1998), jet rotation (e.g., 
ICercmeira fc de Gouveia Dal Pinol|2004 ). jet inlet precession 
" Caproni fc Abraham! l2004al lbT) 



them (e.g.. lHardee et al.ll200lD . 



or a combination of 



Concerning jet precession, lAbraham fc Carrara! (Il998l ) 
developed an analytical model to study precession 
in the quasar 3C 279, which w as also applied to 
the q uasars 3C 273 and 3C345 l|Abraham fc Romero! 
1 19991 : ICaproni fc Abraham! 12004a! ) . the BL Lac object 
OJ287 and the Seyfert g alaxy 3C 120 (|Abrahaml |2000| ; 
I Caproni fc Abraham 2004b). The precession model param- 
eters in those works were estimated from the position an- 
gles, velocities, and epoch of formation of the superluminal 
features in the radio jet, as well as from the long-term peri- 
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odic variability at optical wavelengths in the case of OJ 287, 
3C 345 and 3C 120. 

In order to study jet precession, it is necessary to mon- 
itor the motion of jet knots, collecting data during a long 
enough period of time (two or more precession periods), with 
the interval between consecutive observations as short as 
possible. The reason for this is to avoid misidentification of 
the jet components, which can lead to the wrong determi- 
nation of their kinematic parameters (e.g., proper motion), 
as well as to the wrong estimation of the precession param- 
eters. However, note that such ideal monitoring is hardly 
achieved in practical situations. 

Aiming to overcome, or at least minimise the referred 
problems related to the jet precession analyses, we introduce 
a powerful statistical technique recently developed to deal 
with multi-extremal problems involving optimization: the 
cross-entropy method (hereafter CE). 

The CE analysis was originally used in the optimiza- 
tion of complex co mputer simulation models involving rare 
events simulations j[Rubinstemlll997n . having been modified 
bv lRubinsteinl (|l999h to deal with continuous multi-extremal 
and discrete combinatorial optimization problems. Its the- 
oretical asymptotic conver gence has been demonstrated 



bv iMargolinl l|2004h . while Kroese. Porotskv fc Rubinstein! 
(2006) studied the efficiency of the CE method in solving 
continuous multi-extremal optimization problems. Some ex- 
amples of ro bustness of th e CE m ethod in several situations 
are listed in lde Boer et "a] (|2005h . 

The basic procedures involved in the CE op- 
timization can be summarised as follows (e.g., 
iKroese. Porotskv fc Rubinsteidl2006h : 

(i) Random generation of the initial parameter sample, 
obeying predefined criteria; 

(ii) Selection of the best samples based on some mathe- 
matical criterion; 

(iii) Random generation of updated parameter samples 
from the previous best candidates to be evaluated in the 
next iteration; 

(iv) Optimization process repeats steps (ii) and (iii) until 
a pre-specified stopping criterion is fulfilled. 

In this work we validate our CE jet precession model 
from a variety of benchmark tests built from synthetic pro- 
cessing jet components. They show the great capability of 
our method to determine jet kinematic parameters without 
any identification scheme for the jet knots, as usually made 
in the literature. This paper is structured as follows: in § 
2, we introduce the CE algorithm and the jet precession 
model used in the calculations. Validation tests and the CE 
parameters are also presented in this section. The optimiza- 
tion results for each benchmark test of § 2, as well as for 
three complementary tests are discussed in § 3. Conclusions 
are presented in § 4. 



2 THE CROSS-ENTROPY METHOD AND 
THE JET PRECESSION 

In this section we introduce the cross-entropy method, de- 
scribing how to estimate the precession model parameters 
from the right ascension and declination offsets between the 




Figure 1. Jet and counter jet precession in the source's reference 
frame (see text for the meaning of the symbols). The axis 
points to the observer. 



jet components and the unresolved core obtained from the 
model fittings of radio interferometric observations. 



2.1 Cross-entropy algorithm for continuous 
optimization 

Let us suppose that we wish to study a set of Nd observa- 
tional data in terms of an analytical model characterized by 
N p parameters pi,p 2 ,..., Pn p ■ 

The main goal of the CE continuous multi-extremal 
optimization method is to find the set of parame- 
ters x* = (p*,P2, ■■■iPNp) f° r which the model pro- 
vides the best description of the d ata (jRubins tcin 1999; 
IKroese. Porotskv fc Rubinsteinll2006l ). It is performed gen- 
crating randomly N independent sets of model parameters 
X = (xi,x 2 , xjv), where Xi = (pu,P2i, —,PN P i), and min- 
imizing an objective function S(x) used to transmit the qual- 
ity of the fit during the run process. If the convergence to 
the exact solution is achieved then S(x*) — > 0. 

In order to find the optimal solution from CE optimiza- 
tion, we start by defining the parameter range in which the 
algorithm will search for the best candidates: p™ 1 " ^ Pj(k) ^ 
p™ ax , where k represents the iteration number. Introducing 



Pj(0) = ( P r n +Pr X )/2 and ,7,(0) 
can compute X(0) from: 

X I ,(0)=p J (0)+a J (0)G lJ , 



(pj» _ pf™)/2, we 



(1) 



where Gij is an N x N p matrix with random numbers gen- 
erated from a zero-mean normal distribution with standard 
deviation of unity. 

The next step is to calculate Si(0) for each set of Xj(0), 
ordering them according to increasing values of Si. Then 
the first iVciitJ3 se t of parameters is selected, i.e. the A^iite- 
samples with lowest ,5-values, which will be labeled as the 
elite sample array X elltG (0). 

We then determine the mean and standard deviation of 
the elite sample, pf ltB (0) and cr| llte (0) respectively, as: 



1 The estimation of the optimal or near-optimal value for the 
parameter AT e ij te for our validation tests is discussed in Section 3. 
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Table 1. Precession model parameters used to generate ar- 
tificially the right-ascension and declination offsets of the jet 
components. The sense of precession was chosen to be clock- 
wise, with the precession period at the observer's reference 
frame corresponding to 14.20 yr. 



1 C 



Test 



no 

(deg) 



<f>0 
(deg) 



90 
(deg) 



I ^noise I 
(mas) 



Tl 


8.99 


-161.3 


11.4 


8.7 


3.954 


0.0 


T2 


8.99 


-161.3 


11.4 


8.7 


3.954 


0.1 


T3 


8.99 


-161.3 


11.4 


8.7 


3.954 


0.5 


Tl 


10.01 


-165.0 


15.3 


5.0 


3.971 


0.0 


T5 


7.47 


118.1 


4.0 


0.0 


3.860 


0.0 



—elite / r\\ 

Pj (o) = 



^elite 

E 

i=l 



3 (0), 



2 (o) 



1=1 



(2) 



(3) 



The array X at the next iteration is determined as: 
X ij {l)=0f(O) + <rf a '{O)G ij , (4) 

This process is repeated from equation (2) , with Gij 
regenerated at each iteration. The optimization stops when 
either the mean value of af lte (k) is smaller than a predefined 
value or the maximum number of iterations fc max is reached. 

In order to prevent convergence to a sub-optimal 
solution due to the intrinsic rapid converg e nce o f the 
CE method, iKroese. Porotskv fc RubinstehT |2006) sug- 
gested the implementation of a fixed smoothing scheme for 

— clitc.s/T \ i elite. s / 1 \ 

p i (fc) and a- (k): 



"0 

—elite.; 
elite, s 



(k) 



3 

I -elite 



(fc) 



1 pf te (fc 



n ,"" '"(k) = ctd(k)af lte (k) ■ 



[1 - a d (k)] 



-1) 

elite 



(fc 



(5) 
(6) 



where a' is a smoothing constant parameter (0 < a' < 1) 
and Qd(fc) is a dynamic smoothing parameter at feth itera- 
tion: 



oy(fc) 



(1-k-y 



(7) 



wit h < a < 1 and q is an integer typica lly between 5 and 
10 (|Kroese. Porotskv fc Rubinsteid[2006l ). 

As mentioned before, such parametrisation prevents the 
algorithm from finding a non-global minimum solution since 
it guarantee s polynomial speed of convergence i nstead of 
exponential l|Kroese. Porotskv fc R ubinstein l2006h - 



2.2 Ballistic jet precession model 

Let us consider a relativistic jet receding from the core with 
a constant bulk velocity /3 (in units of the light speed c) 
that precesses around a fixed axis z s , as shown in Fig. 1. 
Jet and counter-jet recede from the nucleus in opposite di- 
rections. Because of precession, the jet inlet direction varies 
with time with a precession period P pr oc, s measured in the 
source reference frame, producing a cone with semi-aperture 
angle ipo. The precession phase co s At s = 27r(t s — io,s)/-P P rec,s 
is chosen arbitrarily to be zero in the j/ s z s -plane at t s — to,s- 
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Figure 2. Jet parameter evolution in the observer's reference 
frame for the precession model tests Tl, T2 and T3. Open squares 
represent the kinematic parameters used to produce the synthetic 
components. 



From these definitions, the instantaneous jet inlet direction 
is given in terms of a unit vector with Cartesian components: 

e x ,s(i s ) = sin^o sin[o; s (t s - £o,s)], 

e y ,s(t s ) = sin tpo cos[io s {t B — to, s )}, 

e z ,s(ts) = cosy?o- 

Making two consecutive clockwise rotations of the 
source's coordinate system, the first one by an angle 0o 
around the y s -axis, and the second by an angle rjo of 
the resulting system arou nd the z-axis, we have (e.g., 
ICaproni fc AbraharrJl2004al lbh: 



e x ,obs(i s ) = A(t B ) cos 770 - e y , s (t s ) sin rjo, 
e y ,obs(£s) = A(t 3 ) sin 770 + e y , s (i s ) cos r/o, 



(8) 
(9) 
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Figure 3. Jet parameter evolution in the observer's reference 
frame for the precession model test T4. Open squares represent 
the kinematic parameters used to produce the synthetic compo- 
nents. 
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Figure 4. Jet parameter evolution in the observer's reference 
frame for the precession model test T5. Open squares represent 
the kinematic parameters used to produce the synthetic compo- 
nents. 



e z ,obs(£ s ) = — e x , s (t s ) sin <f>o + e z , B (t s ) cos cj>o, 
where: 

A(t s ) = e x , 8 (ts) cos 0o + e z , B (t s ) sin <f> . 



(10) 



(11) 



The parameter <f>o is the angle between the precession 
cone axis (positive z s -axis direction) and the line of sight 
(positive direction), and rjo is the position angle of 

the axis on the plane of the sky (positive from north to east 
direction). 

The instantaneous angle between the jet and the line of 
sight </>(t s ) is calculated from: 



while the position angle of the jet on the plane of the sky is 
obtained from: 



v(t. 



arctan 



,obs (^s 



s ) = arccos[e z ,obs(£ B )L 



(12) 



Cx,obs (is 

The observed jet velocity /3 bs(i s ) is: 
/3 bs(i s ) = 7/3<5 (ts) sin 4>(t 3 ), 
where the jet Lorentz factor 7 is 

7=(1-/3T 1/2 , 

and the jet Doppler boosting factor 8 is 
S(t e )=^- 1 [1-/3 cos (fits)}' 1 - 



(13) 

(14) 
(15) 
(16) 



Cross-entropy optimiser 5 



o 

-2 
-4 
-6 



2 0-2 
Aa (mas) 




-3.0 



-2.5 



-2.0 



-1.5 



1.0 -0.5 0.0 



Figure 5. Right ascension and declination offsets, as well as the core-component distances of the synthetic jet components used in 
Tl-test. 



In order to compare predictions from the precession 
model with the observational data, it is necessary to trans- 
form the elapsed time measured in the source's reference 
frame dt s to th e time interval in t he observer's framework 
dt bs- Based on lGower et al.l(jl982l ). we can write the trans- 
formation as: 



At ob 



(17) 



-Pprcc.obs <5 _1 (T)GfT 

where At h B = (iobs — io.obs) and At s = At s /P prcCjS . 

The relation between the precession period in the 
source's framework and that measured by the observer 



prec,obs is given as: 

q *y -fprec,obs 



prcc,£ 



(18) 



where z is the redshift of the source. 

2.3 Building synthetic jet components 

In order to validate our CE precession model, we constructed 
five sets of data, simulating right ascension and declination 
offsets between core and jet components, Aa mo[ i and A<5 mo d 
respectively. These artificial jet knots were generated from 
different sets of precession model parameters, as listed in Ta- 
ble 1, assuming that they recede ballistically from the core 
after their formation. We also imposed arbitrarily a maxi- 
mum angular separation of about 10 mas to the components. 
The evolution of cj>, /3 bs, T], and S in the observer's reference 
frame, as well as /3 bs vs. r\ for the precession models pro- 
vided in Table 1 are presented in Figs. 2-4. 

Tests T1-T3 consist of a relativistic jet whose precession 
produces large and rapid time variations in its position angle 
on the plane of the sky and in its apparent velocity. Changes 
in these quantities in test T4 are smoother, even though the 
jet is more relativistic. Test T5 represents the case of a mild 
relativistic jet with no precession (<^o = 0). 

In the case of Tl, we generated 15 data points that 
trace the trajectories of seven components ejected at the 
epochs shown by the squares in Fig. 2. These data points 
are sparsely sampled in an interval of about 1.5 precession 
periods, in order to mimic observations in which the source 
is not well-sampled in time, as can be seen in Fig. 5. 

Assuming that the jet was observed at nine different 
epochs (see Fig. 5), we calculated Aa mo( j and AS mo d from: 



Aa mo d(t bs) = Ar m od(tobs)sin[r7mod(tej)], 
A<5 mod (iobs) = Ar mod (t bs) cos[»7 m od(tej)], 
where: 

A?* m od (iobs ) — | ^obs tej | /^mod (^ej ) , 



(19) 
(20) 

(21) 



T]mod and /imod are respectively the position angle and the 
component's proper motion at ejection time t e j- 

Tests T2 and T3 are similar to Tl except for the intro- 
duction of uniformly-distributed random noise in the right 
ascension and declination offsets, with maximum amplitude 
A nc> ise shown in Table 1. The noise was added to the original 
coordinates as: 



Aa = Aa mo d + A noise , 
AS = A5 mod + A n ? se) 



(22) 
(23) 



where the superscripts RA and DEC were used to distin- 
guish the random noise introduced into right ascension and 
declination respectively. They simulate errors introduced 
into observational data from the model fitting uncertainties 
and frequency-dependent shifts i n core-component distances 
produced by opacity effects (e.g.. felandford fc KonigH Il979l : 
lLobanovll 19981 : iKovalev et al.ll2008l ). 

The usual procedure before modeling jet precession is 
to identify the jet knots and follow them in time to obtain 
their pro per motions and position angles on the plane of the 
sky (e.g., Abraham & Romero 1999; K ellermann et al ] |2004l : 
ICaproni fc Abr aham 2004a b). The lack of close enough suc- 
cessive observations or high rates of jet component produc- 
tion might contribute to an ambiguous component identifi- 
cation, as well as its accurate kinematic behaviour. 

Motivated by this, we generate randomly 50 compo- 
nents spread out in an interval of about 1.3 precession pe- 
riods, as shown in Fig. 6. Each component is unique (one 
data point for each jet component) and obeys kinetically the 
T4 precession model. Note that the plot of core-component 
distance as a function of time in Fig. 6 seems to suggest 
(wrongly!) the existence of about seven components ejected 
with low velocities and receding non-ballistically from the 
core. 

Test T5 consists of a non-precessing jet, in which the 
components move at a position angle of 118? 1 with constant 
velocity (f3 ohB ~ 6.1). We created 13 data points, correspond- 
ing to seven jet components launched at the epochs repre- 
sented in Fig. 4. Plots of right ascension and declination 
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Figure 6. Right ascension and declination offsets, as well as the core-component distances of the synthetic jet components used in 
T4-tcst. 
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Figure 7. Right ascension and declination offsets, as well as the core-component distances of the synthetic jet components used in 
T5-tcst. 



offsets, as well as those for the core-component distances for 
T5, can be seen in Fig. 7. 

2.4 Estimating precession parameters from 
cross-entropy method 

All validation tests were performed with a fixed P pr oc,obs, 
reducing to five the free model parameters to be optimized 
by the CE method (7, rjo, 4>o, 'Po and Ar s |3- The reason for 
this is that precession period is usually constrained by other 
studies, suc h as periodic variatio n s in the continuum spec- 
trum (e.g.. iGiacconi et al .1 1 1 9731 ; iLehto fc Valtonenl 1 19961 ; 
ICaproni fc Abraham 2004a b) and in velocity or/and inten- 
sity of emission lines (e.g |Margonlll984l ; IStirling et aTll2002i ; 
IStorchi-Bergmann et al.ll2003l ). 

In order to obtain the optimal precession parameters 
through the CE optimization method, the values of N, a, 
iV e iite and q must be provided. We decided to use q = 5 for 
all tests in Table 1 in order to guarantee a polynomial speed 
of convergence, as mentioned in Section 2 (see section 3.2 
for some consequences of departures from that value). For 
the other parameters, we have used 10 ^ N 100, 0.6 ^ 
a < 0.9 and 5 ^ N e u te ^ 10. We also adopted «' = 1, which 



2 The inclusion of P prC c,obs as an additional parameter to be op- 
timised brings the necessity to parallelize the algorithm in order 
to avoid very long processing time (the IDL version of the code 
used in this work spends about 20 hours to generate a partic- 
ular entry in Table 2 in a 2.4 GHZ Intel(R) Core(TM) 2 Quad 
processor). This will be implemented in future. 



makes equation (5) not critical to this work (pf ltc ' B (k) = 



-elite 



(k)). 



At each iteration N sets of random precession parame- 
ters are generated, based on equations (4) and (5), which are 
used to determine /3 bs, ?7obs, 5obs and At b s from equations 
(12)-(17). Then Aa mo( j, A<5 mo( j and Ar mo d are calculated 
using equations (19)-(21) and compared quantitatively with 
the data through the objective function. The elite sample, 
characterized by the AgHte-lowest values of S, is chosen to 
be used as input into equations (4) and (5) to produce the 
next set of precession parameters. 

The objective function S(k) transmits to the CE algo- 
rithm the best tentative solutions at the feth iteration. We 
have chosen S(k) as: 

S(k) = T(fc) +J2 {[ScM? + [SsM? + [Sr z (k)} 2 } . (24) 
1=1 

where: 

S ai (fc) = Aa, - Aa modi (fc) , (25) 

S Si (k)=A8i-AS modl (k), (26) 

5 , r ,(fc) = Ar 1 -Ar modi (fc), (27) 

with Ar 2 = Aa 2 +A<5 2 . It is important to emphasise that the 
terms S ai and Sa i strongly constrain the instantaneous jet 
direction in the optimization process, while inclusion of S ri 
provides additional constraint to the modelling of the core- 
component distances besides improving the convergence per- 
formance of the method. The penalty function T(fc) is given 
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Table 2. Results for the jet precession model parameters from cross-entropy method for q = 5 
and fc m ax = 500. The columns 5-9 show the relative errors among the real precession model 
parameters and their three-run mean values obtained from the CE optimisation (subscripts refer 
to the corresponding precession parameters). 
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(%) 


(%) 


(%) 


(%) 








Tl 


10 


5 


0.6 


67.79 


19.76 


288.71 


335.25 


0.1330 


5.161 


± 


0.190 


Tl 


20 


5 


0.(5 


33.59 


9.80 


144.29 


168.02 


0.1139 


2.327 


± 


0.032 


Tl 


30 


5 


0.(5 


14.19 


21.53 


315.91 


359.77 


0.4984 


4.926 


± 


0.127 


Tl 


40 


5 


0.(5 


7.09 


1.24 


6.74 


11.65 


0.0044 


1.151 


± 


0.080 


Tl 


50 


5 


0.(5 


4.81 


0.24 


0.42 


3.20 


0.0020 


0.963 


± 


0.091 


Tl 


70 


5 


0.(5 


1.79 


0.07 


0.05 


1.27 


0.0089 


0.209 


± 


0.000 


Tl 


100 


5 


0.(5 


7.32 


0.26 


0.41 


3.75 


0.0046 


0.265 


± 


0.011 


Tl 


50 


5 


0.7 


1.00 


0.00 


0.36 


0.26 


0.0011 


0.943 


± 


0.091 


Tl 


50 


5 


0.9 


1.45 


0.01 


0.34 


0.68 


0.0019 


0.900 


± 


0.098 


Tl 


50 


10 
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[ 10 if | *<S max (fe) -tpeak \> e 

where i pea k is the epoch at which there is an outburst in the 
light curve, which we associate with the radiation boosting 
of the underlying jet (see e.g. ICaproni fc Abraham! [2004bl 
for more details), tc* max is the epoch at which the Doppler 
boosting factor 8 reaches £ times its maximum value (0 ^ 
£ ^ 1) and e is a constant interpreted as the uncertainty in 
the epoch of occurrence of the maximum of the jet radiation 
boosting in the light curve. In this work, we adopted £ = 0.9 
and e = 4% of the precession period. 

We can see that the choice of S is based on the minimi- 
sation of the quadratic distances between the observational 
data and those produced by the precession model, while 
T(fc) provides an extra constraint in the Ar s parameter. 
It is important to emphasise that the CE method allows the 
inclusion on the calculations of as many penalty functions 
as necessary to improve the convergence of the method. 

The search for the optimal parameters was halted at 
iteration fc max = 500 for all tests listed in Table 1. Good es- 
timates of the precession parameters were obtained at fc max 



in most cases (see next section). An additional motivation 
to limit the number of iterations was checking the conver- 
gence's behaviour of the CE method as a and A^nte have 
their values varied. 

Each validation test was performed three times, having 
the CE values for 7, r]o,(po, fo and Ar s been calculated aver- 
aging their fc max -values weighted by the inverse of S 2 (fc max ), 
where S is averaged over Aeiite-best solutions at fc max . Their 
respective uncertainties were obtained from the calculation 
of their standard deviations. 

Note that all validation tests presented in this work do 
not use any information concerning counter-jet data. In the 
case of micro-quasars and Herbig-Haro sources, the inclu- 
sion of counter-jet data must increase the performance of 
the CE optimisation. The reason for this is that the pres- 
ence of the counter-jet in the images imposes additional con- 
straint on the Doppler boosting parameter, which on the 
other hand depends on the jet speed and the jet viewing 
angle. On the other hand, the relativistic counter-jets in BL 
Lac objects and quasars are significantly Doppler-deboosted 
for small viewing angles, which difficult their direct detec- 
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tions, as well as any kinematic study of their knots. For radio 
galaxies, the long timescales involved on the changes of the 
jet/counter-jet structures limit the detectability of the pre- 
cession in terms of the component motions. Nevertheless, if 
the jet-counter-jet intensity ratio is known, it can be used 
in all cases above to put additional constraints on the value 
of Doppler-boosting factor during the optimisation process. 



of a provides a faster convergence since the mix between 
consecutive iterations is diminished and the parameters al- 
ways converged to the real model values, while N c y ltc — 10 
increased the convergence time. The best performance was 
achieved for AT e Hte = 5 and a — 0.9 for the tests shown 
in Fig. 9; however we adopted the more conservative value 
a = 0.6 in the remaining tests. 



3 RESULTS AND DISCUSSION 



3.1 Tests T1-T5 



There is not yet a theoretical foundation that guaran- 
tees which CE parameters must be used in the calcula- 
tions in order to maximise the optim isation efficiency (e.g., 
iKroese. Porotskv fe Rubinsteinl2006T ) . Therefore, it is neces- 
sary to explore the CE parameter-space in order to estimate 
the optimal or near-optimal values. 

We present in Table 2 the relative difference between 
the precession model parameters and those derived from the 
CE technique for all tests listed in Table 1, as well as the 
CE parameters N, ATeUte, and a used in the calculations. All 
runs were arbitrarily halted at iteration fc max = 500, which 
means that not necessarily the optimisations reached their 
best value in all cases. Nevertheless, the results presented 
here are representative of the great potential of the cross- 
entropy continuous optimization technique as a quantitative 
tool for jet precession inference. 



3.1.1 Tl 

The first seven runs of Tl were made fixing iVeiite = 5 and 
a = 0.6 varying N from 10 to 100. The precession parame- 
ters were allowed to vary in a wide range : 2.3 ^ 7 ^ 22.4, 
-200° < 770 < -100°, 0?1 < 4> < 50°, 0° < tpo < 40° and 
3.8 ^ At s ^ 4.0. The results are shown in Fig. 8, in which 
the evolution of the precession parameters and the objective 
function are plotted as a function of the iteration number. 

We can see that the increase of N results in a better 
convergence to the correct precession parameters. Indeed, 
only runs with N > 50 were successful in recovering the 
model parameters within 3<r (or 5% in terms of relative er- 
ror). For N < 50, the capability of finding the right solution 
was substantially depreciated due to either the large varia- 
tion of the parameters within the three-run optimization or 
the rapid convergence to sub-optimal parameters (e.g. see 
the evolution of 4>o for iV = 30 in Fig. 



In a general context, Ide Boer et all l|2005h and 
IKroese. Porotskv fc Rubinstein! ( 20061 ) adopted respectively 



N = (5 — 10)JVp and N = 100Ap in their calculations, which 
in our case corresponds to 125 ^ N ^ 250 and N = 500, 
respectively. Our results suggest the existence of a minimum 
value for the sample size N for which the CE optimization 
can be applied to find the jet precession model parameter^]: 
N > WN P . 

We carried out five additional runs for Tl fixing N — 50 
and assuming N e me = 5 and 10, as well as a = 0.6, 0.7 and 
0.9, as shown in Table 2. The performance of the CE op- 
timization is presented in Fig. 9. As expected, the increase 



It is valid for all validation tests presented in this work. 



3.1.2 T2 

T2 validation tests are similar to those of Tl, except for the 
introduction of random noise with maximum amplitude of 
0.1 mas. The objective was to analyse the CE optimization 
performance in a more realistic situation. 

First, we checked the influence of A^ in the optimization 
process, fixing JVeuta = 5 and a = 0.6 and assuming N = 50, 
70 and 100. The results are shown in Table 2. Although 
N — 70 has minimised slightly better the mean value of S, 
the estimates of the precession parameters for A^ = 100 are 
closer to the real ones (< 3%). 

We also carried out five additional runs for T2 fixing 
N — 50 and assuming A^eiite = 5 and 10, as well as a = 0.6, 
0.7 and 0.9. The results are shown in Table 2. Analogously 
to Tl, the best performance has been achieved for AT e ute — 5 
and a = 0.9. 



3.1.3 T3 

Here we added to Tl data random noise with maximum 
amplitude of 0.5 mas. This is an extreme validation test 
since uncertainties in right ascension and declination offsets 
are often smaller, especially at high-frequency observations 
(> 22 GHz), and at core-component distances smaller than 
about 3 mas. The CE optimization behaviour for T3, using 
A'ditc = 5 and a — 0.6 and assuming N — 50, 70 and 100 is 
presented in Table 2. 

The CE technique for N = 50 recovered the precession 
parameters with less efficiency than for TV = 70 and 100, 
reaching relative errors of about 71% in cj>o, 61% in tpo and 
32% in 7. There were no substantial differences in the pre- 
cession parameters values obtained from A^ = 70 and 100 
(< 0.3%), with the relative error smaller than ~ 27% in the 
worst case. 



3.1.4 T4 

Test T4 was thought to be representative of those obser- 
vations in which the jet components are not well-identified 
to be followed along the different epochs of observation. As 
discussed previously, each one of the 50 points in T4 corre- 
sponds to a different jet component. 

This challenging test was conducted maintaining 
A^eiite = 5 and a = 0.6 fixed and assuming N — 50, 70 and 
100 in order to study the influence of N in the optimiza- 
tion process. The results are listed in Table 2. We can see 
that the precession parameters were recovered in all cases, 
with a relative error smaller than 1%. We believe that one 
of the reasons behind this impressive result is the substan- 
tial increase of the number of the jet components used to 
constrain the CE optimization in contrast to previous tests 
(50 against 15 data points in T1-T3). 
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Figure 8. Evolution of the precession parameters (first five columns from left to right) as well as the value of the objective function (last 
column) for Tl (see Table 1) as a function of the iteration number k, assuming W e i; te = 5, a = 0.6 and q = 5. From top to bottom we 
have N = 10, 20, 30, 40, 50, 70 and 100 respectively (small labels were put into 5-plots to guide the reader). Thick solid lines show the real 
precession model parameters. The mean parameter values (weighted by the squared inverse values of the objective function) recovered 
by the cross-entropy algorithm after three independent runs, as well as their 3<r-deviations, are displayed by thin solid and dotted lines 
respectively. 



A similar good performance was obtained for the five 
additional runs for T4, now fixing iV = 50 and assuming 
JVciitc = 5 and 10, as well as a = 0.6, 0.7 and 0.9. The dif- 
ferences in the final precession parameters from those tests 
were no larger 2%, having a slightly better performance been 
observed for iVeiite = 5 and a = 0.7, as it can be seen in Ta- 
ble 2. 

It is important to emphasise that such impressive per- 
formance was obtained using data points without any de- 
viation from the pure predictions from the T4 precession 
model. In practice, jet components could present more com- 
plex motions due to e.g. jet-environment interaction, hydro- 



and magnetohydrodynamical instabilities, as well as core- 
component distance shifts introduced by opacity effects, 
which would decrease the recovering accuracy of the CE 
precession method. Indeed, we made an extra test, adding 
random noise with maximum amplitude of 0.1 mas in the 
right ascension and declination offsets in the 50 data points 
of test T4, assuming N = 100, AT e i ite = 5 and a — 0.6. As 
expected, there was a substantial reduction in the algorithm 
performance, which recovered the real precession parameters 
within a mean relative error of about 13%. However, more 
realistic errors (e.g. smaller for shorter distances to the core) 
are expected to increase the algorithm performance. 
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Figure 9. Evolution of the precession parameters (first five columns from left to right) as well as the value of the objective function (last 
column) for Tl (see Table 1) as a function of the iteration number k, assuming N = 50 and q = 5. From top to bottom, the first set of 
graphs correspond to iV e i; te = 5 and a = 0.6, 0.7 and 0.9 respectively, while for the last set iV e i ite = 10 and a = 0.6, 0.7 and 0.9 (small 
labels were put into 5-plots to guide the reader). Thick solid lines show the real precession model parameters. The mean parameter 
values (weighed by the squared inverse values of the objective function) recovered by the cross-entropy algorithm after three independent 
runs, as well as their 3<r-deviations, are displayed by thin solid and dotted lines respectively. 



3.1.5 T5 

What would be the behaviour of the CE algorithm if there 
were no jet precession? Would the algorithm be able to pro- 
vide any information about this non-precessing jet? To ad- 
dress these questions we performed test T5, a non-precessing 
mildly relativistic jet model with 7 = 7.5 (see Table 1 for 
more details). 

The results from this test are presented in Table 2. The 
parameters 770 and ipo were recovered by the optimization 
process strikingly well. The relative errors in 7 were about 
38%, 68% and 12% for N = 50,70 and 100 respectively, 
while for Ar s they were always smaller than 1%. However the 
CE optimization was not able to find the correct value of 4>q. 
It is important to emphasise that we do not expect to recover 
all parameters in this situation since we are applying a jet 
precession model to a non-precessing jet, which obviously 
can bring some difficulty to the algorithm. 



We believe that the failure in recovering the param- 
eters 4>o and 7 was caused by degenerate solutions when 
the jet viewing angle is constant in time. If <po = 0, /j mo d 
and ?7 m od in equations (19)-(21) are time independent, so 
that for a given value of t e j it is possible to find pi mo d (that 
depends on 7 and 4>o) that provide an acceptable solution 
in the right ascension-declination plane. For example, if we 
decrease |t bs — t c - s \ by a factor two and increase ^i mo d by 
the same factor, the displacement of the component remains 
completely unchanged. An alternative to break this degen- 
eracy would be to introduce an additional constraint on 
Doppler boosting factor (e.g., using some value for Doppler 
factor derived from high-energy observations) , which means 
further constraints to 7 and 0o- Another possibility could be 
the usage of the CE method to verify whether there is any 
signature of jet precession (from whether or not ipo = 0). If 
not, then we could return to the standard procedures based 
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Table 3. Jet precession model parameters from cross-entropy 
method for T2-test considering additional values of q. The CE 
parameters used in the calculations were N = 50, N e n te = 5 and 

^max — 500. 
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69.3 


18.4 
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332.0 
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15.91 ± 0.56 
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16.07 ± 0.90 
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64.1 


16.5 
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289.4 
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9.58 ± 0.14 


10 
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68.6 


17.0 


242.7 


283.9 


0.3 


2.71 ± 0.04 


10 


0.7 


51.5 


12.1 


166.0 


191.0 


0.3 


4.41 ± 0.21 


10 


0.9 


64.9 


16.9 


256.6 


294.4 


0.6 


2.77 ± 0.04 



on some identification component scheme, such as the time 
evolution of the jet knots in a plot of core-component dis- 
tance, which constrains their apparent proper motions and 
hence the values of 7 and 4>Q- 

Nevertheless, it is important to emphasise that the CE 
precession method was capable to find easily the right zero- 
value for ipo, which strongly suggests its usage to verify the 
existence of any precession signature in the observational 
data. 

3.2 Additional tests 

We now discuss additional validation tests concerning 
changes in the values of the CE parameter q and of the 
precession period in the observer's reference frame P P roc,obs- 
The influence of the sense of jet precession on the optimiza- 
tion process will be also discussed in this section. 

3.2.1 Dynamical smoothing 

All validation tests discussed in last section were done as- 
suming q — 5. Roughly speaking, this parameter is responsi- 
ble for co ntrolling the speed of convergence of t he CE opti- 
mization l|Kroese. Porotskv fc Rubinsteinll2006h . Increasing 
its value makes convergence faster, even though it might lead 
to a sub-optimal solution. On the other hand, low values of 
q imply in a slower convergence, requiring higher number of 
iterations to obtain good solutions. 

The precession parameters for additional runs using q = 
1 and 10 with the T2 data are listed in Table 3. For each 
of them, we also varied a to probe possible differences in 
the optimization performance. As expected, the adoption of 
q = 1 and 10 worsened substantially the performance of the 
CE method, so that none of the precession parameters were 
correctly recovered. Furthermore, we can note that S(k lnii ^) 
listed in Table 3 are higher than those obtained in the same 
circumstances assuming q = 5 (see Table 2). 

Therefore, it seems that q « 5 makes the CE optimiza- 
tion more stable and efficient at least in the case of our jet 
precession model. 

3.2.2 Sense of precession 

Tests T1-T4 were constructed adopting clockwise precession 
in the models and in the optimisation procedure. The ques- 
tion that arises is what would happened if we had chosen 
the wrong sense of precession? 



Table 4. Jet precession model parameters from cross-entropy 
method for T2-test considering additional values for the preces- 
sion period (/ is given in units of Pp rec ,obs used to build tests TI- 
TS). The CE parameters used in the calculations were N = 50, 
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8.19 ± 0.37 
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0.34 


8.93 


7.19 


0.1675 


1.674 ± 0.075 
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119.0 


32.9 


13.37 ± 0.04 


2.0 


148.5 


0.4 


154.8 


15.9 


64.4 


14.16 ± 0.00 



In order to provide some quantitative answer to this 
question we used the T2-test data considering the opposite 
direction for jet precession (uj s < 0). The CE parameters 
used in the calculations were N = 50, N e n tc = 5, a = 0.6, 
q — 5 and fc max = 500, which led to 7 = 2.69 ± 0.03, 
770 = -179.83 ± 1.28, 0o = 20.60 ± 1.90, tp = 14.12 ± 1.39, 
Ar s = 3.9288 ± 0.0033 and S(fc max ) = 6.907 ± 0.294. The 
counterclockwise values are far from the real ones and they 
also produce a higher S(fc max ) in comparison to the clock- 
wise precession. Therefore, our CE precession model can also 
indicate the correct sense of the jet precession, if included 
as a CE parameter. 



3.2.3 Precession period 

All optimization tests performed until now used a fixed pe- 
riod P proc ,obs, the same used to build our artificial jet com- 
ponents. But what would happen if we had chosen another 
value for P pr cc,obs? 

To tackle this question, we made four optimization runs 
for T2 data with precession periods corresponding to 0.5, 
0.7,1.4, and 2.0 of the original value. The CE adopted pa- 
rameters were N = 50, N e n ta = 5, a = 0.6 and q = 5. The 
optimization results for the four periods are presented in 
Table 4. 

We can see that the minimum value of S^fcmax) is ob- 
tained for the true period; therefore, in a real situation in 
which this parameter is not well-constrained, we can run 
the algorithm for different values of the precession period, 
choosing the solutions that minimises .S^fcmax). Note that the 
inclusion of the precession period as an additional parameter 
to be optimised by the algorithm lengthens the completion of 
the calculations. Our guess is that it might increase the du- 
ration of each run by at least a factor of 10. This very rough 
estimate is based on a set of CE parameters as adopted in 
this work. 



3.3 Kinematics of the jet components 

In Figs. 10-19 we show the model predictions for the core- 
component distances, and the right ascension and declina- 
tion offsets for the T1-T5 tests. Solid lines in these graphs 
represent the predictions from the CE optimisation adopt- 
ing N — 100, Aeiite = 5, a = 0.6 and q — 5. Concerning 
core-component distance plots, the zero-separation epochs 
refer to those marked by the squares in Figs. 2-4, while the 
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Figure 10. Core-component distances for test Tl. Diamonds rep- 
resent the synthetic jet components, while lines are the predicted 
proper motions for each component found from CE precession 
model optimization. The ejection epochs correspond exactly to 
those shown in Fig. 2 
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Figure 11. Snapshots of right ascension and declination offsets 
of the jet components in test Tl (observation epochs are shown 
on the top of the plots). Diamonds represent the synthetic jet 
components, while lines are the predicted instantaneous jet ori- 
entation provided by the CE precession model optimization. 



inclinations of the lines correspond to the proper motions ob- 
tained from the optimisation processes. In the case of right 
ascension-declination plots, the continuous lines represent 
the snapshot of the precession helices generated from the 
derived precession model parameters. 

Note that these calculations did not use any predefined 
identificati on scheme for the jet knots, as done in previous 
works (e.g., I Abraham fc Carrarall 19981 ; ICaproni fc Abraham! 
l2004al lbh. In other words, it was not provided to the al- 
gorithm any information about which pairs of right ascen- 
sion and declination offsets belong to a given component. It 
means that misidentification of the jet components due to 
the lack of a good time sampling of the observations, which 
might lead to the wrong determination of their kinematic pa- 
rameters, is completely avoided in our method. Therefore, 
it relaxes moderately the necessity of collecting data during 
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Figure 12. The same as Fig. 10 but for the precession test T2. 
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Figure 13. The same as Fig. 11 but for the precession test T2. 



a long period of time, with the interval between consecutive 
observations as short as possible. 

A good example of how misleading may be the iden- 
tification of the jet components only based on kinematic 
features is test T4 presented in Fig. 6. Due to the sparse 
time coverage, we are induced (specially from the core- 
component distance plot!) to assume the existence of 7 to 10 
non-ballistic jet knots, depending on the data association we 
made. Obviously their kinetic parameters would differ com- 
pletely from those calculated considering the correct num- 
ber of (ballistic) components. It is important to emphasise 
that we are not claiming the nonexistence of non-ballistic 
jet motions which have been reported in the literature (e.g . , 
iHoman et all 1200 ll ; Ijorstad et a"i1l2005l ; lAgudo et ai]|2007l ). 
However, some fraction of the bent motions might be the 
result of a misinterpretation of the observational data due 
mainly to their bad sampling in time, as the test T4 sug- 
gests. 
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Figure 14. The same as Fig. 10 but for the precession test T3. 
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Figure 15. The same as Fig. 11 but for the precession test T3. 



4 CONCLUSIONS 

Many astrophysical jets present signatures of precession, 
which were inferred from jet kinematics an d/or continuum 
and lin e-profile periodic va riabili ty (e.g.. | Giacconi et al. 
19731; iLehto fc Valtonenl Il99fj ; ICaproni fe Abraham 
2004al lbT). In the case of jet kinematics, it is necessary 
to determine previously the apparent proper motions of 
the knots, as well as their trajectories on the plane of 
the sky that, on the other hand, depend on the correct 
identification of the jet components. 

To deal with those potential difficulties, we have 
implemented the cross-entropy method for continuous 
multi-extremal optimization with dynamic smoothing 



(K rocsc, Porotskv fc Ru binstein 20061) into o ur jet preces- 
sion model (e.g., Caproni fc Abrahamll2004al ). As far as we 



know, it is the first case of application of the CE technique 
in astrophysics. 

In this method, tentative sets of precession model pa- 
rameters are generated randomly from a normal distribution 
in each iteration, being the best candidates (elite sample) 
used to construct the next set of solutions. The elite sample 
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Figure 16. Core-component distances for test T4. Diamonds rep- 
resent the synthetic jet components, while lines are the predicted 
proper motions for each component found from CE precession 
model optimization. The ejection epochs correspond exactly to 
those shown in Fig. 3 
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Figure 17. The same as Fig. 11 but for the precession test T4. 



was selected from the A'dito-set of parameters that better 
minimises an objective function S. 

In order to validate our CE precession model, we con- 
structed five sets of data, simulating right ascension and 
declination offsets between core and jet components, Acv mo d 
and A5 mo d respectively. These artificial jet components were 
generated from different combinations of precession model 
parameters, assuming that they recede ballistically from the 
core after their ejection. For each of these tests, the CE pre- 
cession method was applied to the data, aiming to recover 
the precession parameters used to create them. 
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Figure 18. Core-component distances for test T5. Diamonds rep- 
resent the synthetic jet components, while lines are the predicted 
proper motions for each component found from CE precession 
model optimization. The ejection epochs correspond exactly to 
those shown in Fig. 4 
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Figure 19. The same as Fig. 11 but for the precession test T5. 



Our results have shown that even in the most challeng- 
ing tests, the CE method was able to find the correct pa- 
rameters within 1%-level in most cases. The models included 
highly relativistic precessing jets, as well as a mildly rel- 
ativistic jet with no precession. The ability to recover the 
correct precession model parameters was dependent on the 
cross-entropy parameters used in the optimization: all vali- 
dation tests performed in this work seems to indicate that 
N > 50, a > 0.7 and q ~ 5 improve the algorithm's ef- 
ficiency. It is important to emphasise that such mapping 
is necessary since there is no theoretical proposition that 
points surely to which CE parameters must be used in the 
calculations to maximise the optimis ation's efficiency (e.g., 
iKroese. Porotskv fe Rubinsteidl2006h . 

Component identification problems related to fitting 
procedures, as well as observations sparsely sampled in time 
and sources with high rate production of jet components, 
may influence the follow up of the jet knots, which conse- 
quently might contribute to a kinematic misinterpretation 



of the data. In order to verify how sensitive is our algorithm 
to this issue, we performed the T4-test that consists of 50 
precessing jet components generated to mimic eight non- 
ballistic knots. Our technique was able to recover the origi- 
nal jet precession parameters, ignoring any previous compo- 
nents misidentification. Note that the inclusion of random 
fluctuations in the right ascension and declination offsets 
decreases the ability of the algorithm in recovering the real 
precession parameters. 

To keep the numerical execution of the algorithm within 
acceptable intervals, we fixed the precession period in the ob- 
server's reference frame and the sense of precession (clock- 
wise or counterclockwise rotation) during the optimization 
processes. This practical limitation led us to perform addi- 
tional tests to check the behaviour of the CE method when 
a wrong choice for the precession period or the sense of ro- 
tation is made. As expected, such additional runs clearly 
indicated that the minimum value of the objective function 
is obtained when the correct values for the precession period 
and sense of rotation are employed. 

Our method was successful in pointing out the lack of 
precession in the case of test T5 (mildly relativistic non- 
precessing jet) , recovering the values of tpo and rjo with a very 
impressive precision. Although reasonable values have been 
also obtained for 7 and Ar s , the method failed completely 
to find the correct value for <j>o. This is not supposed to be a 
serious issue since we are modelling a non-precessing jet in 
terms of a precession model, which obviously can produce 
some inconsistencies in the results. Nevertheless, we argue 
that our CE method can be useful in indicating the presence 
of any precession signature in the observational data. 

Our results show the great potentiality of the CE 
method to study precession (or the lack of) in astrophysical 
jets, utilising only the time variation of the right ascension 
and declination coordinates of the knots. 
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